<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html>
 <head>
  <title>strsna.f</title>
 <meta name="generator" content="emacs 21.3.1; htmlfontify 0.20">
<style type="text/css"><!-- 
body { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.default   { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.default a { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: underline; }
span.string   { color: rgb(188, 143, 143);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.string a { color: rgb(188, 143, 143);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: underline; }
span.comment   { color: rgb(178, 34, 34);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.comment a { color: rgb(178, 34, 34);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: underline; }
 --></style>

 </head>
  <body>

<pre>
      SUBROUTINE <a name="STRSNA.1"></a><a href="strsna.f.html#STRSNA.1">STRSNA</a>( JOB, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
     $                   LDVR, S, SEP, MM, M, WORK, LDWORK, IWORK,
     $                   INFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  -- LAPACK routine (version 3.1) --
</span><span class="comment">*</span><span class="comment">     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
</span><span class="comment">*</span><span class="comment">     November 2006
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Modified to call <a name="SLACN2.9"></a><a href="slacn2.f.html#SLACN2.1">SLACN2</a> in place of <a name="SLACON.9"></a><a href="slacon.f.html#SLACON.1">SLACON</a>, 7 Feb 03, SJH.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     .. Scalar Arguments ..
</span>      CHARACTER          HOWMNY, JOB
      INTEGER            INFO, LDT, LDVL, LDVR, LDWORK, M, MM, N
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Array Arguments ..
</span>      LOGICAL            SELECT( * )
      INTEGER            IWORK( * )
      REAL               S( * ), SEP( * ), T( LDT, * ), VL( LDVL, * ),
     $                   VR( LDVR, * ), WORK( LDWORK, * )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Purpose
</span><span class="comment">*</span><span class="comment">  =======
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  <a name="STRSNA.25"></a><a href="strsna.f.html#STRSNA.1">STRSNA</a> estimates reciprocal condition numbers for specified
</span><span class="comment">*</span><span class="comment">  eigenvalues and/or right eigenvectors of a real upper
</span><span class="comment">*</span><span class="comment">  quasi-triangular matrix T (or of any matrix Q*T*Q**T with Q
</span><span class="comment">*</span><span class="comment">  orthogonal).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  T must be in Schur canonical form (as returned by <a name="SHSEQR.30"></a><a href="shseqr.f.html#SHSEQR.1">SHSEQR</a>), that is,
</span><span class="comment">*</span><span class="comment">  block upper triangular with 1-by-1 and 2-by-2 diagonal blocks; each
</span><span class="comment">*</span><span class="comment">  2-by-2 diagonal block has its diagonal elements equal and its
</span><span class="comment">*</span><span class="comment">  off-diagonal elements of opposite sign.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Arguments
</span><span class="comment">*</span><span class="comment">  =========
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  JOB     (input) CHARACTER*1
</span><span class="comment">*</span><span class="comment">          Specifies whether condition numbers are required for
</span><span class="comment">*</span><span class="comment">          eigenvalues (S) or eigenvectors (SEP):
</span><span class="comment">*</span><span class="comment">          = 'E': for eigenvalues only (S);
</span><span class="comment">*</span><span class="comment">          = 'V': for eigenvectors only (SEP);
</span><span class="comment">*</span><span class="comment">          = 'B': for both eigenvalues and eigenvectors (S and SEP).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  HOWMNY  (input) CHARACTER*1
</span><span class="comment">*</span><span class="comment">          = 'A': compute condition numbers for all eigenpairs;
</span><span class="comment">*</span><span class="comment">          = 'S': compute condition numbers for selected eigenpairs
</span><span class="comment">*</span><span class="comment">                 specified by the array SELECT.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  SELECT  (input) LOGICAL array, dimension (N)
</span><span class="comment">*</span><span class="comment">          If HOWMNY = 'S', SELECT specifies the eigenpairs for which
</span><span class="comment">*</span><span class="comment">          condition numbers are required. To select condition numbers
</span><span class="comment">*</span><span class="comment">          for the eigenpair corresponding to a real eigenvalue w(j),
</span><span class="comment">*</span><span class="comment">          SELECT(j) must be set to .TRUE.. To select condition numbers
</span><span class="comment">*</span><span class="comment">          corresponding to a complex conjugate pair of eigenvalues w(j)
</span><span class="comment">*</span><span class="comment">          and w(j+1), either SELECT(j) or SELECT(j+1) or both, must be
</span><span class="comment">*</span><span class="comment">          set to .TRUE..
</span><span class="comment">*</span><span class="comment">          If HOWMNY = 'A', SELECT is not referenced.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  N       (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The order of the matrix T. N &gt;= 0.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  T       (input) REAL array, dimension (LDT,N)
</span><span class="comment">*</span><span class="comment">          The upper quasi-triangular matrix T, in Schur canonical form.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  LDT     (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The leading dimension of the array T. LDT &gt;= max(1,N).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  VL      (input) REAL array, dimension (LDVL,M)
</span><span class="comment">*</span><span class="comment">          If JOB = 'E' or 'B', VL must contain left eigenvectors of T
</span><span class="comment">*</span><span class="comment">          (or of any Q*T*Q**T with Q orthogonal), corresponding to the
</span><span class="comment">*</span><span class="comment">          eigenpairs specified by HOWMNY and SELECT. The eigenvectors
</span><span class="comment">*</span><span class="comment">          must be stored in consecutive columns of VL, as returned by
</span><span class="comment">*</span><span class="comment">          <a name="SHSEIN.74"></a><a href="shsein.f.html#SHSEIN.1">SHSEIN</a> or <a name="STREVC.74"></a><a href="strevc.f.html#STREVC.1">STREVC</a>.
</span><span class="comment">*</span><span class="comment">          If JOB = 'V', VL is not referenced.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  LDVL    (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The leading dimension of the array VL.
</span><span class="comment">*</span><span class="comment">          LDVL &gt;= 1; and if JOB = 'E' or 'B', LDVL &gt;= N.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  VR      (input) REAL array, dimension (LDVR,M)
</span><span class="comment">*</span><span class="comment">          If JOB = 'E' or 'B', VR must contain right eigenvectors of T
</span><span class="comment">*</span><span class="comment">          (or of any Q*T*Q**T with Q orthogonal), corresponding to the
</span><span class="comment">*</span><span class="comment">          eigenpairs specified by HOWMNY and SELECT. The eigenvectors
</span><span class="comment">*</span><span class="comment">          must be stored in consecutive columns of VR, as returned by
</span><span class="comment">*</span><span class="comment">          <a name="SHSEIN.86"></a><a href="shsein.f.html#SHSEIN.1">SHSEIN</a> or <a name="STREVC.86"></a><a href="strevc.f.html#STREVC.1">STREVC</a>.
</span><span class="comment">*</span><span class="comment">          If JOB = 'V', VR is not referenced.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  LDVR    (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The leading dimension of the array VR.
</span><span class="comment">*</span><span class="comment">          LDVR &gt;= 1; and if JOB = 'E' or 'B', LDVR &gt;= N.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  S       (output) REAL array, dimension (MM)
</span><span class="comment">*</span><span class="comment">          If JOB = 'E' or 'B', the reciprocal condition numbers of the
</span><span class="comment">*</span><span class="comment">          selected eigenvalues, stored in consecutive elements of the
</span><span class="comment">*</span><span class="comment">          array. For a complex conjugate pair of eigenvalues two
</span><span class="comment">*</span><span class="comment">          consecutive elements of S are set to the same value. Thus
</span><span class="comment">*</span><span class="comment">          S(j), SEP(j), and the j-th columns of VL and VR all
</span><span class="comment">*</span><span class="comment">          correspond to the same eigenpair (but not in general the
</span><span class="comment">*</span><span class="comment">          j-th eigenpair, unless all eigenpairs are selected).
</span><span class="comment">*</span><span class="comment">          If JOB = 'V', S is not referenced.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  SEP     (output) REAL array, dimension (MM)
</span><span class="comment">*</span><span class="comment">          If JOB = 'V' or 'B', the estimated reciprocal condition
</span><span class="comment">*</span><span class="comment">          numbers of the selected eigenvectors, stored in consecutive
</span><span class="comment">*</span><span class="comment">          elements of the array. For a complex eigenvector two
</span><span class="comment">*</span><span class="comment">          consecutive elements of SEP are set to the same value. If
</span><span class="comment">*</span><span class="comment">          the eigenvalues cannot be reordered to compute SEP(j), SEP(j)
</span><span class="comment">*</span><span class="comment">          is set to 0; this can only occur when the true value would be
</span><span class="comment">*</span><span class="comment">          very small anyway.
</span><span class="comment">*</span><span class="comment">          If JOB = 'E', SEP is not referenced.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  MM      (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The number of elements in the arrays S (if JOB = 'E' or 'B')
</span><span class="comment">*</span><span class="comment">           and/or SEP (if JOB = 'V' or 'B'). MM &gt;= M.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  M       (output) INTEGER
</span><span class="comment">*</span><span class="comment">          The number of elements of the arrays S and/or SEP actually
</span><span class="comment">*</span><span class="comment">          used to store the estimated condition numbers.
</span><span class="comment">*</span><span class="comment">          If HOWMNY = 'A', M is set to N.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  WORK    (workspace) REAL array, dimension (LDWORK,N+6)
</span><span class="comment">*</span><span class="comment">          If JOB = 'E', WORK is not referenced.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  LDWORK  (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The leading dimension of the array WORK.
</span><span class="comment">*</span><span class="comment">          LDWORK &gt;= 1; and if JOB = 'V' or 'B', LDWORK &gt;= N.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  IWORK   (workspace) INTEGER array, dimension (2*(N-1))
</span><span class="comment">*</span><span class="comment">          If JOB = 'E', IWORK is not referenced.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  INFO    (output) INTEGER
</span><span class="comment">*</span><span class="comment">          = 0: successful exit
</span><span class="comment">*</span><span class="comment">          &lt; 0: if INFO = -i, the i-th argument had an illegal value
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Further Details
</span><span class="comment">*</span><span class="comment">  ===============
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  The reciprocal of the condition number of an eigenvalue lambda is
</span><span class="comment">*</span><span class="comment">  defined as
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">          S(lambda) = |v'*u| / (norm(u)*norm(v))
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  where u and v are the right and left eigenvectors of T corresponding
</span><span class="comment">*</span><span class="comment">  to lambda; v' denotes the conjugate-transpose of v, and norm(u)
</span><span class="comment">*</span><span class="comment">  denotes the Euclidean norm. These reciprocal condition numbers always
</span><span class="comment">*</span><span class="comment">  lie between zero (very badly conditioned) and one (very well
</span><span class="comment">*</span><span class="comment">  conditioned). If n = 1, S(lambda) is defined to be 1.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  An approximate error bound for a computed eigenvalue W(i) is given by
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                      EPS * norm(T) / S(i)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  where EPS is the machine precision.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  The reciprocal of the condition number of the right eigenvector u
</span><span class="comment">*</span><span class="comment">  corresponding to lambda is defined as follows. Suppose
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              T = ( lambda  c  )
</span><span class="comment">*</span><span class="comment">                  (   0    T22 )
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Then the reciprocal condition number is
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">          SEP( lambda, T22 ) = sigma-min( T22 - lambda*I )
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  where sigma-min denotes the smallest singular value. We approximate
</span><span class="comment">*</span><span class="comment">  the smallest singular value by the reciprocal of an estimate of the
</span><span class="comment">*</span><span class="comment">  one-norm of the inverse of T22 - lambda*I. If n = 1, SEP(1) is
</span><span class="comment">*</span><span class="comment">  defined to be abs(T(1,1)).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  An approximate error bound for a computed right eigenvector VR(i)
</span><span class="comment">*</span><span class="comment">  is given by
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                      EPS * norm(T) / SEP(i)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  =====================================================================
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     .. Parameters ..
</span>      REAL               ZERO, ONE, TWO
      PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0, TWO = 2.0E+0 )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Scalars ..
</span>      LOGICAL            PAIR, SOMCON, WANTBH, WANTS, WANTSP
      INTEGER            I, IERR, IFST, ILST, J, K, KASE, KS, N2, NN
      REAL               BIGNUM, COND, CS, DELTA, DUMM, EPS, EST, LNRM,
     $                   MU, PROD, PROD1, PROD2, RNRM, SCALE, SMLNUM, SN
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Arrays ..
</span>      INTEGER            ISAVE( 3 )
      REAL               DUMMY( 1 )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Functions ..
</span>      LOGICAL            <a name="LSAME.193"></a><a href="lsame.f.html#LSAME.1">LSAME</a>
      REAL               SDOT, <a name="SLAMCH.194"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>, <a name="SLAPY2.194"></a><a href="slapy2.f.html#SLAPY2.1">SLAPY2</a>, SNRM2
      EXTERNAL           <a name="LSAME.195"></a><a href="lsame.f.html#LSAME.1">LSAME</a>, SDOT, <a name="SLAMCH.195"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>, <a name="SLAPY2.195"></a><a href="slapy2.f.html#SLAPY2.1">SLAPY2</a>, SNRM2
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Subroutines ..
</span>      EXTERNAL           <a name="SLABAD.198"></a><a href="slabad.f.html#SLABAD.1">SLABAD</a>, <a name="SLACN2.198"></a><a href="slacn2.f.html#SLACN2.1">SLACN2</a>, <a name="SLACPY.198"></a><a href="slacpy.f.html#SLACPY.1">SLACPY</a>, <a name="SLAQTR.198"></a><a href="slaqtr.f.html#SLAQTR.1">SLAQTR</a>, <a name="STREXC.198"></a><a href="strexc.f.html#STREXC.1">STREXC</a>, <a name="XERBLA.198"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Intrinsic Functions ..
</span>      INTRINSIC          ABS, MAX, SQRT
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Executable Statements ..
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Decode and test the input parameters
</span><span class="comment">*</span><span class="comment">
</span>      WANTBH = <a name="LSAME.207"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOB, <span class="string">'B'</span> )
      WANTS = <a name="LSAME.208"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOB, <span class="string">'E'</span> ) .OR. WANTBH
      WANTSP = <a name="LSAME.209"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOB, <span class="string">'V'</span> ) .OR. WANTBH
<span class="comment">*</span><span class="comment">
</span>      SOMCON = <a name="LSAME.211"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( HOWMNY, <span class="string">'S'</span> )
<span class="comment">*</span><span class="comment">
</span>      INFO = 0
      IF( .NOT.WANTS .AND. .NOT.WANTSP ) THEN
         INFO = -1
      ELSE IF( .NOT.<a name="LSAME.216"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( HOWMNY, <span class="string">'A'</span> ) .AND. .NOT.SOMCON ) THEN
         INFO = -2
      ELSE IF( N.LT.0 ) THEN
         INFO = -4
      ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
         INFO = -6
      ELSE IF( LDVL.LT.1 .OR. ( WANTS .AND. LDVL.LT.N ) ) THEN
         INFO = -8
      ELSE IF( LDVR.LT.1 .OR. ( WANTS .AND. LDVR.LT.N ) ) THEN
         INFO = -10
      ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Set M to the number of eigenpairs for which condition numbers
</span><span class="comment">*</span><span class="comment">        are required, and test MM.
</span><span class="comment">*</span><span class="comment">
</span>         IF( SOMCON ) THEN
            M = 0
            PAIR = .FALSE.
            DO 10 K = 1, N
               IF( PAIR ) THEN
                  PAIR = .FALSE.
               ELSE
                  IF( K.LT.N ) THEN
                     IF( T( K+1, K ).EQ.ZERO ) THEN
                        IF( SELECT( K ) )
     $                     M = M + 1
                     ELSE
                        PAIR = .TRUE.
                        IF( SELECT( K ) .OR. SELECT( K+1 ) )
     $                     M = M + 2
                     END IF
                  ELSE
                     IF( SELECT( N ) )
     $                  M = M + 1
                  END IF
               END IF
   10       CONTINUE
         ELSE
            M = N
         END IF
<span class="comment">*</span><span class="comment">
</span>         IF( MM.LT.M ) THEN
            INFO = -13
         ELSE IF( LDWORK.LT.1 .OR. ( WANTSP .AND. LDWORK.LT.N ) ) THEN
            INFO = -16
         END IF
      END IF
      IF( INFO.NE.0 ) THEN
         CALL <a name="XERBLA.264"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="STRSNA.264"></a><a href="strsna.f.html#STRSNA.1">STRSNA</a>'</span>, -INFO )
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Quick return if possible
</span><span class="comment">*</span><span class="comment">
</span>      IF( N.EQ.0 )
     $   RETURN
<span class="comment">*</span><span class="comment">
</span>      IF( N.EQ.1 ) THEN
         IF( SOMCON ) THEN
            IF( .NOT.SELECT( 1 ) )
     $         RETURN
         END IF
         IF( WANTS )
     $      S( 1 ) = ONE
         IF( WANTSP )
     $      SEP( 1 ) = ABS( T( 1, 1 ) )
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Get machine constants
</span><span class="comment">*</span><span class="comment">
</span>      EPS = <a name="SLAMCH.287"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>( <span class="string">'P'</span> )
      SMLNUM = <a name="SLAMCH.288"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>( <span class="string">'S'</span> ) / EPS
      BIGNUM = ONE / SMLNUM
      CALL <a name="SLABAD.290"></a><a href="slabad.f.html#SLABAD.1">SLABAD</a>( SMLNUM, BIGNUM )
<span class="comment">*</span><span class="comment">
</span>      KS = 0
      PAIR = .FALSE.
      DO 60 K = 1, N
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Determine whether T(k,k) begins a 1-by-1 or 2-by-2 block.
</span><span class="comment">*</span><span class="comment">
</span>         IF( PAIR ) THEN
            PAIR = .FALSE.
            GO TO 60
         ELSE
            IF( K.LT.N )
     $         PAIR = T( K+1, K ).NE.ZERO
         END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Determine whether condition numbers are required for the k-th
</span><span class="comment">*</span><span class="comment">        eigenpair.
</span><span class="comment">*</span><span class="comment">
</span>         IF( SOMCON ) THEN
            IF( PAIR ) THEN
               IF( .NOT.SELECT( K ) .AND. .NOT.SELECT( K+1 ) )
     $            GO TO 60
            ELSE
               IF( .NOT.SELECT( K ) )
     $            GO TO 60
            END IF
         END IF
<span class="comment">*</span><span class="comment">
</span>         KS = KS + 1
<span class="comment">*</span><span class="comment">
</span>         IF( WANTS ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Compute the reciprocal condition number of the k-th
</span><span class="comment">*</span><span class="comment">           eigenvalue.
</span><span class="comment">*</span><span class="comment">
</span>            IF( .NOT.PAIR ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Real eigenvalue.
</span><span class="comment">*</span><span class="comment">
</span>               PROD = SDOT( N, VR( 1, KS ), 1, VL( 1, KS ), 1 )
               RNRM = SNRM2( N, VR( 1, KS ), 1 )
               LNRM = SNRM2( N, VL( 1, KS ), 1 )
               S( KS ) = ABS( PROD ) / ( RNRM*LNRM )
            ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Complex eigenvalue.
</span><span class="comment">*</span><span class="comment">
</span>               PROD1 = SDOT( N, VR( 1, KS ), 1, VL( 1, KS ), 1 )
               PROD1 = PROD1 + SDOT( N, VR( 1, KS+1 ), 1, VL( 1, KS+1 ),
     $                 1 )
               PROD2 = SDOT( N, VL( 1, KS ), 1, VR( 1, KS+1 ), 1 )
               PROD2 = PROD2 - SDOT( N, VL( 1, KS+1 ), 1, VR( 1, KS ),
     $                 1 )
               RNRM = <a name="SLAPY2.344"></a><a href="slapy2.f.html#SLAPY2.1">SLAPY2</a>( SNRM2( N, VR( 1, KS ), 1 ),
     $                SNRM2( N, VR( 1, KS+1 ), 1 ) )
               LNRM = <a name="SLAPY2.346"></a><a href="slapy2.f.html#SLAPY2.1">SLAPY2</a>( SNRM2( N, VL( 1, KS ), 1 ),
     $                SNRM2( N, VL( 1, KS+1 ), 1 ) )
               COND = <a name="SLAPY2.348"></a><a href="slapy2.f.html#SLAPY2.1">SLAPY2</a>( PROD1, PROD2 ) / ( RNRM*LNRM )
               S( KS ) = COND
               S( KS+1 ) = COND
            END IF
         END IF
<span class="comment">*</span><span class="comment">
</span>         IF( WANTSP ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Estimate the reciprocal condition number of the k-th
</span><span class="comment">*</span><span class="comment">           eigenvector.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Copy the matrix T to the array WORK and swap the diagonal
</span><span class="comment">*</span><span class="comment">           block beginning at T(k,k) to the (1,1) position.
</span><span class="comment">*</span><span class="comment">
</span>            CALL <a name="SLACPY.362"></a><a href="slacpy.f.html#SLACPY.1">SLACPY</a>( <span class="string">'Full'</span>, N, N, T, LDT, WORK, LDWORK )
            IFST = K
            ILST = 1
            CALL <a name="STREXC.365"></a><a href="strexc.f.html#STREXC.1">STREXC</a>( <span class="string">'No Q'</span>, N, WORK, LDWORK, DUMMY, 1, IFST, ILST,
     $                   WORK( 1, N+1 ), IERR )
<span class="comment">*</span><span class="comment">
</span>            IF( IERR.EQ.1 .OR. IERR.EQ.2 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Could not swap because blocks not well separated
</span><span class="comment">*</span><span class="comment">
</span>               SCALE = ONE
               EST = BIGNUM
            ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Reordering successful
</span><span class="comment">*</span><span class="comment">
</span>               IF( WORK( 2, 1 ).EQ.ZERO ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Form C = T22 - lambda*I in WORK(2:N,2:N).
</span><span class="comment">*</span><span class="comment">
</span>                  DO 20 I = 2, N
                     WORK( I, I ) = WORK( I, I ) - WORK( 1, 1 )
   20             CONTINUE
                  N2 = 1
                  NN = N - 1
               ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Triangularize the 2 by 2 block by unitary
</span><span class="comment">*</span><span class="comment">                 transformation U = [  cs   i*ss ]
</span><span class="comment">*</span><span class="comment">                                    [ i*ss   cs  ].
</span><span class="comment">*</span><span class="comment">                 such that the (1,1) position of WORK is complex
</span><span class="comment">*</span><span class="comment">                 eigenvalue lambda with positive imaginary part. (2,2)
</span><span class="comment">*</span><span class="comment">                 position of WORK is the complex eigenvalue lambda
</span><span class="comment">*</span><span class="comment">                 with negative imaginary  part.
</span><span class="comment">*</span><span class="comment">
</span>                  MU = SQRT( ABS( WORK( 1, 2 ) ) )*
     $                 SQRT( ABS( WORK( 2, 1 ) ) )
                  DELTA = <a name="SLAPY2.399"></a><a href="slapy2.f.html#SLAPY2.1">SLAPY2</a>( MU, WORK( 2, 1 ) )
                  CS = MU / DELTA
                  SN = -WORK( 2, 1 ) / DELTA
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Form
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 C' = WORK(2:N,2:N) + i*[rwork(1) ..... rwork(n-1) ]
</span><span class="comment">*</span><span class="comment">                                        [   mu                     ]
</span><span class="comment">*</span><span class="comment">                                        [         ..               ]
</span><span class="comment">*</span><span class="comment">                                        [             ..           ]
</span><span class="comment">*</span><span class="comment">                                        [                  mu      ]
</span><span class="comment">*</span><span class="comment">                 where C' is conjugate transpose of complex matrix C,
</span><span class="comment">*</span><span class="comment">                 and RWORK is stored starting in the N+1-st column of
</span><span class="comment">*</span><span class="comment">                 WORK.
</span><span class="comment">*</span><span class="comment">
</span>                  DO 30 J = 3, N
                     WORK( 2, J ) = CS*WORK( 2, J )
                     WORK( J, J ) = WORK( J, J ) - WORK( 1, 1 )
   30             CONTINUE
                  WORK( 2, 2 ) = ZERO
<span class="comment">*</span><span class="comment">
</span>                  WORK( 1, N+1 ) = TWO*MU
                  DO 40 I = 2, N - 1
                     WORK( I, N+1 ) = SN*WORK( 1, I+1 )
   40             CONTINUE
                  N2 = 2
                  NN = 2*( N-1 )
               END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Estimate norm(inv(C'))
</span><span class="comment">*</span><span class="comment">
</span>               EST = ZERO
               KASE = 0
   50          CONTINUE
               CALL <a name="SLACN2.433"></a><a href="slacn2.f.html#SLACN2.1">SLACN2</a>( NN, WORK( 1, N+2 ), WORK( 1, N+4 ), IWORK,
     $                      EST, KASE, ISAVE )
               IF( KASE.NE.0 ) THEN
                  IF( KASE.EQ.1 ) THEN
                     IF( N2.EQ.1 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                       Real eigenvalue: solve C'*x = scale*c.
</span><span class="comment">*</span><span class="comment">
</span>                        CALL <a name="SLAQTR.441"></a><a href="slaqtr.f.html#SLAQTR.1">SLAQTR</a>( .TRUE., .TRUE., N-1, WORK( 2, 2 ),
     $                               LDWORK, DUMMY, DUMM, SCALE,
     $                               WORK( 1, N+4 ), WORK( 1, N+6 ),
     $                               IERR )
                     ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                       Complex eigenvalue: solve
</span><span class="comment">*</span><span class="comment">                       C'*(p+iq) = scale*(c+id) in real arithmetic.
</span><span class="comment">*</span><span class="comment">
</span>                        CALL <a name="SLAQTR.450"></a><a href="slaqtr.f.html#SLAQTR.1">SLAQTR</a>( .TRUE., .FALSE., N-1, WORK( 2, 2 ),
     $                               LDWORK, WORK( 1, N+1 ), MU, SCALE,
     $                               WORK( 1, N+4 ), WORK( 1, N+6 ),
     $                               IERR )
                     END IF
                  ELSE
                     IF( N2.EQ.1 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                       Real eigenvalue: solve C*x = scale*c.
</span><span class="comment">*</span><span class="comment">
</span>                        CALL <a name="SLAQTR.460"></a><a href="slaqtr.f.html#SLAQTR.1">SLAQTR</a>( .FALSE., .TRUE., N-1, WORK( 2, 2 ),
     $                               LDWORK, DUMMY, DUMM, SCALE,
     $                               WORK( 1, N+4 ), WORK( 1, N+6 ),
     $                               IERR )
                     ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                       Complex eigenvalue: solve
</span><span class="comment">*</span><span class="comment">                       C*(p+iq) = scale*(c+id) in real arithmetic.
</span><span class="comment">*</span><span class="comment">
</span>                        CALL <a name="SLAQTR.469"></a><a href="slaqtr.f.html#SLAQTR.1">SLAQTR</a>( .FALSE., .FALSE., N-1,
     $                               WORK( 2, 2 ), LDWORK,
     $                               WORK( 1, N+1 ), MU, SCALE,
     $                               WORK( 1, N+4 ), WORK( 1, N+6 ),
     $                               IERR )
<span class="comment">*</span><span class="comment">
</span>                     END IF
                  END IF
<span class="comment">*</span><span class="comment">
</span>                  GO TO 50
               END IF
            END IF
<span class="comment">*</span><span class="comment">
</span>            SEP( KS ) = SCALE / MAX( EST, SMLNUM )
            IF( PAIR )
     $         SEP( KS+1 ) = SEP( KS )
         END IF
<span class="comment">*</span><span class="comment">
</span>         IF( PAIR )
     $      KS = KS + 1
<span class="comment">*</span><span class="comment">
</span>   60 CONTINUE
      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End of <a name="STRSNA.493"></a><a href="strsna.f.html#STRSNA.1">STRSNA</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>
